
uppertime<- 6
lowertime<- -6


eventmeans<-FALSE

source("./psid_eventpreamble.R")

#Empyet must be redefined to allow checking of control units too.
rawdata[,empyet:=NULL]

rawdata[,haspartner := marnow]
rawdata[,marnow:=as.numeric(marit==1)]

rawdata[onset_age_hsev>=onset_age + 5,onset_age:=NA]
rawdata<-rawdata[!is.na(onset_age),]

rawdata[,lag_marnow:=c(NA,marnow[-.N]),by=newid]
rawdata[is.na(lag_marnow),lag_marnow:=marnow]
rawdata[,marchange:=cumsum(lag_marnow!=marnow),by=newid]
rawdata[,marchange:=marchange>=1]

for(var in c("ly","soc_sec","trall_nonss","nonhy","y")){
  rawdata[,eval(var):=get(var)/1000]
}

covariates<-1
timevar = "age"
cohortvar = "onset_age"
anycohortvar = "anyonset_age"




setorder(rawdata,newid,year)
rawdata[,onSSI:=(tr_ssi > 0 & !is.na(tr_ssi))]
rawdata[,id:=NULL]

rawdata[,empyet:=cumsum(emp)>=1,by=newid]	
rawdata[,empyet:=empyet == 1 & sample_analytic == 1,by=newid]
rawdata[,empyet_young:=empyet == 1 & sample_analytic == 1 & age <=35,by=newid]
rawdata[,empyet_lowA := empyet == 1 & w_no_bus_eq <= 10]
rawdata[is.na(w_no_bus_eq),empyet_lowA := 0]
rawdata[,empyet_fixmar:=empyet == 1 & marchange == 0]


rawdata[,currentmar:=marnow]
rawdata[,marnow1:=c(marnow[-1],NA),by=newid]
rawdata[is.na(marnow1),marnow1:=marnow]
rawdata[,marnow:=marnow1]
rawdata[,marnow1:=NULL]

rawdata[,foodonly:=totfood]
rawdata[,foodonly_eq:=totfood_eq]
rawdata[,logfoodonly:=log(foodonly)]
rawdata[,logfoodonly_eq:=log(foodonly_eq)]

rawdata[,whiteblack:=white]
rawdata[race==3,whiteblack:=NA]

rawdata_all<-copy(rawdata)
rawdata_all[,highedu:=as.numeric(educ>=3)]
rawdata<-rawdata[educ<=2,]



#Printing descriptive table:
setorder(rawdata,DS)
rawdata[sample_analytic==1 , mean(trouble,na.rm=TRUE),by=DS]
rawdata[sample_analytic==0 , mean(trouble,na.rm=TRUE),by=DS]
rawdata[sample_analytic==1 , mean(diagnostic,na.rm=TRUE),by=DS]
rawdata[sample_analytic==0 , mean(diagnostic,na.rm=TRUE),by=DS]

tab<-
  TR(c("","(1)","(2)","(3)")) +
  TR(c("","Work Limitation Severity"),cspan=c(1,3)) + 
  midrulep(list(c(2,4))) + 
  TR(c("","None","Moderate","Severe"))+
  midrulep(list(c(2,4)))  +
  TR(c("Panel A: Analytic Sample"),cspan=c(4), surround = "\\textit{%s}")+
  TR(c("Has some daily limitation"))%:%TR(c(rawdata[sample_analytic==1 , mean(trouble,na.rm=TRUE),by=DS]$V1),dec = 2) +
  TR(c("Serious health event"))%:%TR(c(rawdata[sample_analytic==1 , mean(diagnostic,na.rm=TRUE),by=DS]$V1),dec = 2) +
 vspace(5)+
 TR(c("Panel B: Persons Excluded due to Non-Work"),cspan=c(4), surround = "\\textit{%s}")+
 TR(c("Has some daily limitation"))%:%TR(c(rawdata[sample_analytic==0 , mean(trouble,na.rm=TRUE),by=DS]$V1),dec = 2) +
 TR(c("Serious health event"))%:%TR(c(rawdata[sample_analytic==0 , mean(diagnostic,na.rm=TRUE),by=DS]$V1),dec = 2)
  
TS(tab, file="health_measures_table",
   output_path="./PSID_eventstudy/",
   pretty_rules=T,
   header=c('r',rep('c',3)))
setorder(rawdata,"marnow")


stat_row<-function(rawdata,varname,variable,dec){
  coeftable<-feols(as.formula(paste0(variable, "~ as.factor(marnow) - 1")), data = rawdata[sample_analytic==1,],cluster="newid")$coeftable
  diff<-feols(as.formula(paste0(variable, "~ marnow")), data = rawdata[sample_analytic==1,],cluster="newid")$coeftable
  TR(varname)%:%TR(coeftable[,"Estimate"],dec=dec,
                   pvalues=coeftable[,"Pr(>|t|)"])%:%
    TR(diff["marnow","Estimate"],dec=dec,
       pvalues=diff["marnow","Pr(>|t|)"])+
    TR("")%:%TR(coeftable[,"Std. Error"],se=TRUE,dec=dec)%:%
    TR(diff["marnow","Std. Error"],se=TRUE,dec=dec)
    
}
rawdata[,n:=1]

  
eventdata<-create_event_data(rawdata[!is.na(anyDI),],
                             timevar=timevar, unitvar="newid", 
                             cohortvar = cohortvar,
                             anycohortvar = anycohortvar,
                             base_time = "baseyear",
                             lower_event_time = lowertime,
                             upper_event_time = uppertime,
                             base_restrict = "empyet",
                             covariate_base_stratify = "marnow",
                             covariate_base_balance = covariates)


eventdata_nowsev <-create_event_data(rawdata[!is.na(anyDI) & (DS== 2 | age < onset_age),],
                                     timevar=timevar, unitvar="newid",
                                     cohortvar = cohortvar,
                                     anycohortvar = anycohortvar,
                                     base_time = "baseyear",
                                     lower_event_time = lowertime,
                                     upper_event_time = uppertime,
                                     base_restrict = "empyet",
                                     covariate_base_stratify = "marnow",
                                   covariate_base_balance = covariates)

eventdata_notsev <-create_event_data(rawdata[!is.na(anyDI) & (DS!= 2 | age < onset_age),],
                                     timevar=timevar, unitvar="newid",
                                     cohortvar = cohortvar,
                                     anycohortvar = anycohortvar,
                                     base_time = "baseyear",
                                     lower_event_time = lowertime,
                                     upper_event_time = uppertime,
                                     base_restrict = "empyet",
                                     covariate_base_stratify = "marnow",
                                     covariate_base_balance = covariates)


eventdata_fixage<-create_event_data(rawdata[!is.na(anyDI),],
                                    timevar=timevar, unitvar="newid", 
                                    cohortvar = cohortvar,
                                    anycohortvar = anycohortvar,
                                    base_time = "baseyear",
                                    lower_event_time = lowertime,
                                    upper_event_time = uppertime,
                                    base_restrict = "empyet",
                             covariate_base_stratify = "marnow",
                             covariate_base_balance = covariates,
                             stratify_balance_val=1,
                             cohortbin=10)

eventdata_chunk<-copy(eventdata)
eventdata_chunk[as.numeric(as.character(event_time))>0 & as.numeric(as.character(event_time)) <= 1,event_time:="0"]
eventdata_chunk[as.numeric(as.character(event_time))>1 ,event_time:="1"]
eventdata_fixage_chunk<-copy(eventdata_fixage)
eventdata_fixage_chunk[as.numeric(as.character(event_time))>0 & as.numeric(as.character(event_time)) <= 1,event_time:="0"]
eventdata_fixage_chunk[as.numeric(as.character(event_time))>1 ,event_time:="1"]

vars<-c("anyDI","currentmar","emp","emp_part","empspouse","empspouse_part","ly","soc_sec","trall_nonss","nonhy","y","y_posttax","hours","sp_hours",
        "ly_eq","soc_sec_eq","trall_nonss_eq","nonhy_eq","y_eq","y_posttax_eq")
varnames<-c("DI Receipt","Currently Married","Full-Time Work","Part-Time Work", "Partner FT Work", "Partner PT Work",
            "Ref. Work Inc.","Soc. Sec. Inc.","Other Transfers","Work Inc. (Others)","Gross Income","Post-Tax Income", "Hours Worked","Partner Hours Worked")
results<-event_ATTs(eventdata,outcomes=vars)
results_nowDI<-event_ATTs(eventdata_nowDI,outcomes=vars)
results_notDI<-event_ATTs(eventdata_notDI,outcomes=vars)

results_chunk<-event_ATTs(eventdata_chunk,outcomes=vars)
results_fixage<-event_ATTs(eventdata_fixage,outcomes=vars,weights="pweight_stratbal")

#checking effects for people who currently are still work-limited:
results_nowsev<-event_ATTs(eventdata_nowsev,outcomes=vars)
results_notsev<-event_ATTs(eventdata_notsev,outcomes=vars)

householdvars<-c("haspartner", "kids", "adults",  "leavenext", "dienext", "dropnext", "endsample")
householdvarnames<-c("Lives with Partner","Number of Kids","Number of Adults","Leaves Sample Household","Dies","Vanishes from Sample","Sample Data Ends")

results_household<-event_ATTs(eventdata,outcomes=householdvars)

results_household_chunk<-event_ATTs(eventdata_chunk,outcomes=householdvars)
results_household_fixage<-event_ATTs(eventdata_fixage,outcomes=householdvars,weights="pweight_stratbal")
results_household_chunk_fixage<-event_ATTs(eventdata_fixage,outcomes=householdvars,weights="pweight_stratbal")

results_household_nowsev<-event_ATTs(eventdata_nowsev,outcomes=householdvars)
results_household_notsev<-event_ATTs(eventdata_notsev,outcomes=householdvars)



foodvars<-c("foodonly_eq","logfoodonly_eq")
foodvarnames<-c("Food Spending","Log Food Spending")
  
eventdata_food<-create_event_data(rawdata[!is.na(totfood_eq) & totfood_eq> 0,],
                                  timevar=timevar, unitvar="newid", 
                                  cohortvar = cohortvar,
                                  anycohortvar = anycohortvar,
                                  base_time = "baseyear",
                                  lower_event_time = lowertime,
                                  upper_event_time = uppertime,
                                  base_restrict = "empyet",
                                  covariate_base_stratify = "marnow",
                                  covariate_base_balance = covariates)


eventdata_food_fixage<-create_event_data(rawdata[!is.na(totfood_eq) & totfood_eq> 0,],
                                  timevar=timevar, unitvar="newid", 
                                  cohortvar = cohortvar,
                                  anycohortvar = anycohortvar,
                                  base_time = "baseyear",
                                  lower_event_time = lowertime,
                                  upper_event_time = uppertime,
                                  base_restrict = "empyet",
                                  covariate_base_stratify = "marnow",
                                  covariate_base_balance = covariates,
                                  stratify_balance_val=1,
                                  cohortbin=10)

results_food<-event_ATTs(eventdata_food,outcomes=foodvars)
results_food_fixage<-event_ATTs(eventdata_food_fixage,outcomes=foodvars,weights="pweight_stratbal")

#Need to do consumption and wealth separately, since they have lots of missing years:
consvars<-c("spending_eq","logspending_eq","totfood_eq","health_eq",
            "education_eq","util_eq","childcare_eq",
            "transport_eq","rent_eq","homeinsure_eq","spending","logspending")
consvarnames<-c("Total Spending","Log Total Spending","Food","Health Care","Education","Utilities","Child Care","Transportation","Rent","Home Insurance")
for(var in consvars){
  rawdata[!is.na(spending_eq)&is.na(get(var)),eval(var):=0]
  rawdata[is.na(spending_eq)&!is.na(get(var)),eval(var):=NA]
}

eventdata_cons<-create_event_data(rawdata[!is.na(spending_eq) & spending_eq> 0,],
                                  timevar=timevar, unitvar="newid", 
                                  cohortvar = cohortvar,
                                  anycohortvar = anycohortvar,
                                  base_time = "baseyear",
                                  lower_event_time = lowertime,
                                  upper_event_time = uppertime,
                                  base_restrict = "empyet",
                                  covariate_base_stratify = "marnow",
                                  covariate_base_balance = covariates)

eventdata_cons_nowsev <-create_event_data(rawdata[!is.na(spending_eq) & spending_eq> 0 & (DS== 2 | age < onset_age),],
                                          timevar=timevar, unitvar="newid",
                                          cohortvar = cohortvar,
                                          anycohortvar = anycohortvar,
                                          base_time = "baseyear",
                                          lower_event_time = lowertime,
                                          upper_event_time = uppertime,
                                          base_restrict = "empyet",
                                  covariate_base_stratify = "marnow",
                                  covariate_base_balance = covariates)

eventdata_cons_notsev <-create_event_data(rawdata[!is.na(spending_eq) & spending_eq> 0 & (DS!= 2 | age < onset_age),],
                                          timevar=timevar, unitvar="newid",
                                          cohortvar = cohortvar,
                                          anycohortvar = anycohortvar,
                                          base_time = "baseyear",
                                          lower_event_time = lowertime,
                                          upper_event_time = uppertime,
                                     base_restrict = "empyet",
                                     covariate_base_stratify = "marnow",
                                     covariate_base_balance = covariates)


eventdata_cons_fixage<-create_event_data(rawdata[!is.na(spending_eq) & spending_eq> 0,],
                                         timevar=timevar, unitvar="newid", 
                                         cohortvar = cohortvar,
                                         anycohortvar = anycohortvar,
                                         base_time = "baseyear",
                                         lower_event_time = lowertime,
                                         upper_event_time = uppertime,
                                  base_restrict = "empyet",
                                  covariate_base_stratify = "marnow",
                                  covariate_base_balance = covariates,
                                  stratify_balance_val=1)

eventdata_cons_chunk<-copy(eventdata_cons)
eventdata_cons_chunk[as.numeric(as.character(event_time))>0 & as.numeric(as.character(event_time)) <= 1,event_time:="0"]
eventdata_cons_chunk[as.numeric(as.character(event_time))>1 ,event_time:="1"]

eventdata_cons_fixage_chunk<-copy(eventdata_cons_fixage)
eventdata_cons_fixage_chunk[as.numeric(as.character(event_time))>0 & as.numeric(as.character(event_time)) <= 1,event_time:="0"]
eventdata_cons_fixage_chunk[as.numeric(as.character(event_time))>1 ,event_time:="1"]

results_cons<-event_ATTs(eventdata_cons,outcomes=consvars)

#checking effects for people who currently are still work-limited:
results_cons_nowsev<-event_ATTs(eventdata_cons_nowsev,outcomes=consvars)
results_cons_notsev<-event_ATTs(eventdata_cons_notsev,outcomes=consvars)


results_cons_chunk<-event_ATTs(eventdata_cons_chunk,outcomes=consvars)
results_cons_fixage<-event_ATTs(eventdata_cons_fixage,outcomes=consvars,weights="pweight_stratbal")

######################

eventdata_cons[as.numeric(as.character(event_time)) < 0 & treated == 1  & pweight > 0 & !is.na(pweight),mean(spending_eq),by=marnow]
summary(results_cons$pooled$spending_eq)


#And now separately for wealth variables:
wealthvars<-c("logw_no_bus_eq","logw_no_bus","logw_no_h_eq","logw_with_h_eq","w_no_bus_eq","w_no_bus","w_no_h_eq","w_with_h_eq")
wealthvarnames<-c("Wealth", "Wealth","Wealth (no Housing)", "Wealth (with Business)")
eventdata_wealth<-create_event_data(rawdata[!is.na(w_no_bus) & w_no_bus> 0 & year >= 2000,],
                                    timevar=timevar, unitvar="newid", 
                                    cohortvar = cohortvar,
                                    anycohortvar = anycohortvar,
                                    base_time = "baseyear",
                                    lower_event_time = lowertime,
                                    upper_event_time = uppertime,
                                    base_restrict = "empyet",
                                    covariate_base_stratify = "marnow",
                                    covariate_base_balance = covariates)

eventdata_wealth_nowsev <-create_event_data(rawdata[!is.na(w_no_bus) & w_no_bus> 0 & year >= 2000 & (DS== 2 | age < onset_age),],
                                            timevar=timevar, unitvar="newid",
                                            cohortvar = cohortvar,
                                            anycohortvar = anycohortvar,
                                            base_time = "baseyear",
                                            lower_event_time = lowertime,
                                            upper_event_time = uppertime,
                                          base_restrict = "empyet",
                                          covariate_base_stratify = "marnow",
                                          covariate_base_balance = covariates)

eventdata_wealth_notsev <-create_event_data(rawdata[!is.na(w_no_bus) & w_no_bus> 0 & year >= 2000 & (DS!= 2 | age < onset_age),],
                                            timevar=timevar, unitvar="newid",
                                            cohortvar = cohortvar,
                                            anycohortvar = anycohortvar,
                                            base_time = "baseyear",
                                            lower_event_time = lowertime,
                                            upper_event_time = uppertime,
                                          base_restrict = "empyet",
                                          covariate_base_stratify = "marnow",
                                          covariate_base_balance = covariates)



eventdata_wealth_fixage<-create_event_data(rawdata[!is.na(w_no_bus) & w_no_bus> 0 & year >= 2000,],
                                           timevar=timevar, unitvar="newid", 
                                           cohortvar = cohortvar,
                                           anycohortvar = anycohortvar,
                                           base_time = "baseyear",
                                           lower_event_time = lowertime,
                                           upper_event_time = uppertime,
                                    base_restrict = "empyet",
                                    covariate_base_stratify = "marnow",
                                    covariate_base_balance = covariates,
                                    stratify_balance_val=1)

eventdata_wealth_chunk<-copy(eventdata_wealth)
eventdata_wealth_chunk[as.numeric(as.character(event_time))>0 & as.numeric(as.character(event_time)) <= 1,event_time:="0"]
eventdata_wealth_chunk[as.numeric(as.character(event_time))>1 ,event_time:="1"]


eventdata_wealth_fixage_chunk<-copy(eventdata_wealth_fixage)
eventdata_wealth_fixage_chunk[as.numeric(as.character(event_time))>0 & as.numeric(as.character(event_time)) <= 1,event_time:="0"]
eventdata_wealth_fixage_chunk[as.numeric(as.character(event_time))>1 ,event_time:="1"]

results_wealth<-event_ATTs(eventdata_wealth,outcomes=wealthvars)


#checking effects for people who currently are still work-limited:
results_wealth_nowsev<-event_ATTs(eventdata_wealth_nowsev,outcomes=wealthvars)
results_wealth_notsev<-event_ATTs(eventdata_wealth_notsev,outcomes=wealthvars)


results_wealth_chunk<-event_ATTs(eventdata_wealth_chunk,outcomes=wealthvars)
results_wealth_fixage<-event_ATTs(eventdata_wealth_fixage,outcomes=wealthvars,weights="pweight_stratbal")
results_wealth_chunk_fixage<-event_ATTs(eventdata_wealth_fixage_chunk,outcomes=wealthvars,weights="pweight_stratbal")

summary(results_wealth$pooled$logw_no_bus_eq)
summary(results_wealth$pooled$logw_with_h_eq)
eventdata_wealth[as.numeric(as.character(event_time)) < 0 & treated == 1  & pweight > 0 & !is.na(pweight),mean(w_no_bus_eq),by=marnow]


#####################################################################
#####################PLOTTING RESULTS################################
decimals<-2


print_chunks<-function(name, results, results_household, results_cons, results_wealth, eventdata, eventdata_cons, eventdata_wealth, weights,noconswealth,nochunk=TRUE){
panelnames<- c("Panel A: Disability Insurance Receipt and Work","Panel B: Household Income","Panel C: Spending and Wealth","","Panel D: Attrition")
if(noconswealth==FALSE){
  vargroups<-list(c("anyDI","hours","sp_hours","currentmar"), 
                  c("ly_eq","soc_sec_eq","trall_nonss_eq","nonhy_eq","y_eq","y_posttax_eq"#,"y_posttax"
                    ),
#                  c("logfoodonly_eq"),
                  c("logspending_eq"),
                  c("logw_no_bus_eq"),
                  c("leavenext", "dienext", "dropnext")
#                  c("logspending"),
#                  c("logw_no_bus")
)
  
  csv_vargroups<-list(c("anyDI","hours","sp_hours","emp","empspouse"), 
                      c("ly_eq","soc_sec_eq","trall_nonss_eq","nonhy_eq","y_eq","y_posttax_eq","y_posttax"),
#                      c("logfoodonly_eq"),
                      c("logspending_eq","spending_eq"),
                      c("logw_no_bus_eq"),
                      c("leavenext", "dienext", "dropnext")
#                      c("spending"),
#                      c("w_no_bus")
)
}
else{
  vargroups<-list(c("anyDI","hours","sp_hours","haspartner"), 
                  c("ly_eq","soc_sec_eq","trall_nonss_eq","nonhy_eq","y_eq","y_posttax_eq"#,"y_posttax"
                    ))
  
  csv_vargroups<-list(c("anyDI","hours","sp_hours","emp","empspouse"), 
                      c("ly_eq","soc_sec_eq","trall_nonss_eq","nonhy_eq","y_eq","y_posttax_eq","y_posttax"))
  
}
vargroupnames<-list(c("DI Receipt","Ref. Hours Worked", "Partner Hours Worked","Currently Married"),
               c("Ref. Work Inc.","Soc. Sec. Inc.","Other Transfers","Work Inc. (Others)","Gross Income","Post-Tax Income"#, "Post-Tax Income (Unequiv.)"
                 ),
#               c("Log Food Spending"),
               c("Log Total Spending"),
               c("Log Wealth"),
               c("Leaves Sample Household","Dies","Vanishes from Sample")
#               c("Log Spending (Unequiv.)"),
#               c("Log Wealth (Unequiv.)"),
               )



csvtab<-NULL
csvtab_dynamic<-NULL
tabheader<-TR(c("","Pre-Disability Means","Effects of Work-Limiting Disability"),cspan=c(1,2,3)) + 
  midrulep(list(c(2,3),c(4,6))) + 
  TR(c("","(1)","(2)","(3)","(4)","(5)")) +
  TR(c("","Single","Partnered","Single","Partnered","(4) - (3)"))+
  midrulep(list(c(2,3),c(4,6))) 
tab<-tabheader

for(varpos in 1:length(vargroups)){
  if(varpos==1 | varpos == 2 ) {
    table_results <- results
    tabledata <- eventdata[as.numeric(as.character(event_time))<0,]
  }
  # if(varpos==3){
  #   table_results <- results_food
  #   tabledata<- eventdata_food[as.numeric(as.character(event_time))<0,]
  # }
  if(varpos==3 ){
    table_results <- results_cons
    tabledata <- eventdata_cons[as.numeric(as.character(event_time))<0,]
  }
  if(varpos==4 ){
    table_results <- results_wealth
    tabledata <- eventdata_wealth[as.numeric(as.character(event_time))<0,]
  }
  if(varpos==5){
    table_results <- results_household
    tabledata<-eventdata[as.numeric(as.character(event_time))<0,]
    
  }
  
  if(varpos != 4) tabnext<- TR(panelnames[varpos],cspan=c(6), surround = "\\textit{%s}")
  for(var in vargroups[[varpos]]){
    namepos <- which(vargroups[[varpos]]==var)
    pos <- which(gsub("lhs: ", "", names(table_results$dynamic))==var)

    means<-feols(as.formula(paste0(var, "~ as.factor(marnow) -1")),
                 data = tabledata[treated == 1 & !is.na(get(weights)),], 
                 weights= tabledata[treated == 1 & !is.na(get(weights)),get(weights)],
                 cluster="id")
    
    tabnext <- tabnext + 
      TR(vargroupnames[[varpos]][namepos])%:%
      TR(c(means$coefficients["as.factor(marnow)0"] , 
           means$coefficients["as.factor(marnow)1"],
           table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.0"]],
           table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.1"]],
           table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.1"]] - table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.0"]]
      ),
      pvalues=c(ifelse(is.na(means$coeftable["as.factor(marnow)0","Pr(>|t|)"]),1,means$coeftable["as.factor(marnow)0","Pr(>|t|)"]) ,
                ifelse(is.na(means$coeftable["as.factor(marnow)1","Pr(>|t|)"]),1,means$coeftable["as.factor(marnow)1","Pr(>|t|)"]),
                ifelse(is.na(table_results$pooled[[pos]]$coeftable["treated_post_stratify1.0","Pr(>|t|)"]),1,table_results$pooled[[pos]]$coeftable["treated_post_stratify1.0","Pr(>|t|)"]),
                ifelse(is.na(table_results$pooled[[pos]]$coeftable["treated_post_stratify1.1","Pr(>|t|)"]),1,table_results$pooled[[pos]]$coeftable["treated_post_stratify1.1","Pr(>|t|)"]),
                2*pt(abs(table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.1"]] - table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.0"]])/(
                  (table_results$pooled[[pos]]$se[["treated_post_stratify1.1"]]^2+table_results$pooled[[pos]]$se[["treated_post_stratify1.0"]]^2)^0.5),
                  df = table_results$pooled[[pos]]$nobs -table_results$pooled[[pos]]$nparams,
                  lower.tail=FALSE)
      ), dec=decimals
      )+
      TR("")%:%
      TR(c(means$se["as.factor(marnow)0"] , means$se["as.factor(marnow)1"],
           table_results$pooled[[pos]]$se[["treated_post_stratify1.0"]],
           table_results$pooled[[pos]]$se[["treated_post_stratify1.1"]],
           (table_results$pooled[[pos]]$se[["treated_post_stratify1.1"]]^2+table_results$pooled[[pos]]$se[["treated_post_stratify1.0"]]^2)^0.5), dec=decimals, se=TRUE)
      vspace(6) 
  }
  if(varpos != 3) tab<-tab+tabnext
  
  if(varpos != 3) TS(tabheader + tabnext, file=paste0(paste(c(name,"pooled_table_p"),collapse="_"),varpos),
                     output_path="./PSID_eventstudy/",
                     pretty_rules=T,
     header=c('r',rep('c',5)))
  
  for(var in csv_vargroups[[varpos]]){
    pos <- which(gsub("lhs: ", "", names(table_results$dynamic))==var)
    row<-t(c(table_results$pooled[[pos]]$coefficients[[c("treated_post_stratify1.0")]],
           table_results$pooled[[pos]]$se[[c("treated_post_stratify1.0")]]))
    rownames(row) <- paste0(var,"_m0")
    csvtab<-rbind(csvtab,row)
  
    row<-t(c(table_results$pooled[[pos]]$coefficients[[c("treated_post_stratify1.1")]],
           table_results$pooled[[pos]]$se[[c("treated_post_stratify1.1")]]))
    rownames(row) <- paste0(var,"_m1")
    csvtab<-rbind(csvtab, row)
  }
  
  if(nochunk==FALSE){
  for(var in csv_vargroups[[varpos]]){
    pos <- which(gsub("lhs: ", "", names(table_results$dynamic))==var)
    row<-t(c(0,table_results$dynamic[[pos]]$coefficients[[c("treated_event_time_stratify0.0")]],
             table_results$dynamic[[pos]]$se[[c("treated_event_time_stratify0.0")]]))
    rownames(row) <- c(paste0(var,"_m0"))
    csvtab_dynamic<-rbind(csvtab_dynamic,row)
    row<-t(c(1,table_results$dynamic[[pos]]$coefficients[[c("treated_event_time_stratify1.0")]],
             table_results$dynamic[[pos]]$se[[c("treated_event_time_stratify1.0")]]))
    rownames(row) <- c(paste0(var,"_m0"))
    csvtab_dynamic<-rbind(csvtab_dynamic,row)
    
    row<-t(c(0,table_results$dynamic[[pos]]$coefficients[[c("treated_event_time_stratify0.1")]],
             table_results$dynamic[[pos]]$se[[c("treated_event_time_stratify0.1")]]))
    rownames(row) <- c(paste0(var,"_m1"))
    csvtab_dynamic<-rbind(csvtab_dynamic, row)
    row<-t(c(1,table_results$dynamic[[pos]]$coefficients[[c("treated_event_time_stratify0.1")]],
             table_results$dynamic[[pos]]$se[[c("treated_event_time_stratify0.1")]]))
    rownames(row) <- c(paste0(var,"_m1"))
    csvtab_dynamic<-rbind(csvtab_dynamic, row)
  }
  }
}

write.csv(csvtab,file=paste0(paste(c(name,"pooled_table.csv"),collapse="_")), 
          row.names = TRUE, col.names = FALSE)

if(nochunk==FALSE) write.csv(csvtab_dynamic,file=paste0(paste(c(name,"dynamic_table.csv"),collapse="_")), 
          row.names = TRUE, col.names = FALSE)


TS(tab, file=paste0(paste(c(name,"pooled_table"),collapse="_")),
   output_path="./PSID_eventstudy/",
   pretty_rules=T,
   header=c('r',rep('c',5)))

table_results <- results_household
tabledata <- eventdata

tab<-TR(c("","Pre-Disability Means","Effects of Work-Limiting Disability"),cspan=c(1,2,3)) + 
  midrulep(list(c(2,3),c(4,6))) + 
  TR(c("","(1)","(2)","(3)","(4)","(5)")) +
  TR(c("","Single","Partnered","Single","Partnered","(4) - (3)"))+
  midrulep(list(c(2,3),c(4,6)))  

for(var in householdvars){
  namepos <- which(householdvars==var)
  pos <- which(gsub("lhs: ", "", names(results_household_chunk$dynamic))==var)

  means<-feols(as.formula(paste0(var, "~ as.factor(marnow) -1")),
               data = tabledata[treated == 1 & !is.na(get(weights)),], 
               weights= tabledata[treated == 1 & !is.na(get(weights)),get(weights)],
               cluster="id")
    
    
    tab <- tab + 
      TR(householdvarnames[namepos])%:%
      TR(c(means$coefficients["as.factor(marnow)0"] , 
           means$coefficients["as.factor(marnow)1"],
           table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.0"]],
           table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.1"]],
           table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.1"]] - table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.0"]]
      ),
      pvalues=c(ifelse(is.na(means$coeftable["as.factor(marnow)0","Pr(>|t|)"]),1,means$coeftable["as.factor(marnow)0","Pr(>|t|)"]) ,
                ifelse(is.na(means$coeftable["as.factor(marnow)1","Pr(>|t|)"]),1,means$coeftable["as.factor(marnow)1","Pr(>|t|)"]),
                ifelse(is.na(table_results$pooled[[pos]]$coeftable["treated_post_stratify1.0","Pr(>|t|)"]),1,table_results$pooled[[pos]]$coeftable["treated_post_stratify1.0","Pr(>|t|)"]),
                ifelse(is.na(table_results$pooled[[pos]]$coeftable["treated_post_stratify1.1","Pr(>|t|)"]),1,table_results$pooled[[pos]]$coeftable["treated_post_stratify1.1","Pr(>|t|)"]),
                2*pt(abs(table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.1"]] - table_results$pooled[[pos]]$coefficients[["treated_post_stratify1.0"]])/(
                  (table_results$pooled[[pos]]$se[["treated_post_stratify1.1"]]^2+table_results$pooled[[pos]]$se[["treated_post_stratify1.0"]]^2)^0.5),
                  df = table_results$pooled[[pos]]$nobs -table_results$pooled[[pos]]$nparams,
                  lower.tail=FALSE)
      ), dec=decimals
      )+
      TR("")%:%
      TR(c(means$se["as.factor(marnow)0"] , means$se["as.factor(marnow)1"],
           table_results$pooled[[pos]]$se[["treated_post_stratify1.0"]],
           table_results$pooled[[pos]]$se[["treated_post_stratify1.1"]],
           (table_results$pooled[[pos]]$se[["treated_post_stratify1.1"]]^2+table_results$pooled[[pos]]$se[["treated_post_stratify1.0"]]^2)^0.5), dec=decimals, se=TRUE)
    vspace(6) 
}

TS(tab, file=paste0(paste(c(name,"householdcomp_table"),collapse="_")),
   output_path="./PSID_eventstudy/",
   pretty_rules=T,
   header=c('r',rep('c',5)))
}


print_chunks("headsev", results, results_household, results_cons, results_wealth, eventdata, eventdata_cons, eventdata_wealth, "pweight",noconswealth=FALSE)
print_chunks("headsev_nowsev", results_nowsev, results_household_nowsev, results_cons_nowsev, results_wealth_nowsev, eventdata_nowsev, eventdata_cons_nowsev, eventdata_wealth_nowsev, "pweight",noconswealth=FALSE)
print_chunks("headsev_notsev", results_notsev, results_household_notsev, results_cons_notsev, results_wealth_notsev, eventdata_notsev, eventdata_cons_notsev, eventdata_wealth_notsev, "pweight",noconswealth=FALSE)
print_chunks("headsev_agebal", results_fixage,  results_household_fixage, results_cons_fixage, results_wealth_fixage, eventdata_fixage, eventdata_cons_fixage, eventdata_wealth_fixage, "pweight_stratbal",noconswealth=FALSE)


#Plot of pre-period means and pre-tests
panelnames<- c("Panel A: Disability Insurance Receipt and Work","Panel B: Household Income","Panel C: Spending and Wealth","", "Panel D: Attrition")
vargroups<-list(c("anyDI","hours","sp_hours","currentmar"), 
                c("ly_eq","soc_sec_eq","trall_nonss_eq","nonhy_eq","y_eq","y_posttax_eq"),
                #c("logfoodonly_eq"),
                c("logspending_eq"),
                c("logw_no_bus_eq"),
                c("leavenext", "dienext", "dropnext")
)
vargroupnames<-list(c("DI Receipt","Ref. Hours Worked","Partner Hours Worked","Currently Married"),
                    c("Ref. Work Inc.","Soc. Sec. Inc.","Other Transfers","Work Inc. (Others)","Gross Income","Post-Tax Income"),
                  #  c("Log Food Spending"),
                    c("Log Total Spending"),
                    c("Log Wealth"),
                    c("Leaves Sample Household","Dies","Vanishes from Sample")
)
tab<-TR(c("","Pre-Trend"),cspan=c(1,2)) + 
  midrulep(list(c(2,3))) +
  TR(c("","(1)","(2)")) +
  TR(c("","Single","Partnered"))+
  midrulep(list(c(2,3)))  

for(varpos in 1:length(vargroups)){
  if(varpos==1 | varpos == 2) {
    table_results <- results_chunk
    tabledata <- eventdata
  }
  if(varpos==3){
    table_results <- results_cons_chunk
    tabledata <- eventdata_cons
  }
  if(varpos==4){
    table_results <- results_wealth_chunk
    tabledata <- eventdata_wealth
  }
  if(varpos==5){
    table_results <- results_household_chunk
    tabledata <- eventdata
  }
  
  if(varpos != 4) tab<- tab + TR(panelnames[varpos],cspan=c(5), surround = "\\textit{%s}")
  for(var in vargroups[[varpos]]){
  namepos <- which(vargroups[[varpos]]==var)
  pos <- which(gsub("lhs: ", "", names(table_results$dynamic))==var)
  
  
  tab<- tab +
    TR(vargroupnames[[varpos]][namepos])%:%
    TR(c(-table_results$pooled[[pos]]$coefficients[["treated_pre_stratify1.0"]],
         -table_results$pooled[[pos]]$coefficients[["treated_pre_stratify1.1"]]
    ),
    pvalues=c(ifelse(is.na(table_results$pooled[[pos]]$coeftable["treated_pre_stratify1.0","Pr(>|t|)"]),1,table_results$pooled[[pos]]$coeftable["treated_pre_stratify1.0","Pr(>|t|)"]),
              ifelse(is.na(table_results$pooled[[pos]]$coeftable["treated_pre_stratify1.1","Pr(>|t|)"]),1,table_results$pooled[[pos]]$coeftable["treated_pre_stratify1.1","Pr(>|t|)"])
    ), dec=decimals
    )+
    TR("")%:%
    TR(c(table_results$pooled[[pos]]$se[["treated_pre_stratify1.0"]],
         table_results$pooled[[pos]]$se[["treated_pre_stratify1.1"]]), dec=decimals, se=TRUE)
  
  }
}



TS(tab, file=paste0(paste(c("headsev","pretrend_table"),collapse="_")),
   output_path="./PSID_eventstudy/",
   pretty_rules=T,
   header=c('r',rep('c',2)))


